%% simulates firms over the lifcycle
i_exit=zeros(n_sim,T);
n_exit=zeros(n_sim,1);
b_exit=zeros(n_sim,1);
z_exit=zeros(n_sim,1);
t_exit=ones(n_sim,1)*T*dt;
inde_z_hist=zeros(n_sim,T);
inde_b_hist=zeros(n_sim,T);
z_hist=zeros(n_sim,T);
b_hist=zeros(n_sim,T);
l_b_hist=zeros(n_sim,T);
f_hist=zeros(n_sim,T);

for ii=1:n_sim
    bd=rand;
    b0=b_up*(bd<=q_hig)+b_md*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn*(bd>q_hig+q_mid);
    z0=zlow*(bd<=q_hig)+zmid*(bd>q_hig)*(bd<=q_hig+q_mid)+zhig*(bd>q_hig+q_mid);

    t=0;inde_alive=1;t_sim=0;
    while t_sim<T && inde_alive==1
        t=t+dt;t_sim=t_sim+1;

        ell_t       =  ppval(pp_ell,max(0,b0));
        drift_I_t   =  ppval(pp_drift,max(0,b0));
        drift_dlb   = dt*(ell_t/b0+sigma^2/2-rho-drift_I_t);
        drift_dlz   = dt*(drift_I_t-sigma^2/2);
        dw  = randn;
        l_b = log(b0) + drift_dlb-sqrt(dt)*sigma*dw;b0  = exp(l_b);
        l_z = log(z0) + drift_dlz+sqrt(dt)*sigma*dw;z0  = exp(l_z);
        if b0<b_bar
            b_hist(ii,t_sim)  =b0;
            l_b_hist(ii,t_sim)=log(b0);
            z_hist(ii,t_sim)  =z0;
        else
            neg = rand;
            if neg >phi || b0*a>b_bar
                inde_alive= 0;
                t_exit(ii)= t;
                i_exit(ii,t_sim)= 1;
            else
                b0=b0*a;
                b_hist(ii,t_sim)  =b0;
                l_b_hist(ii,t_sim)=log(b0);
            end
        end
    end
end

%% yearly averages
tot_years       = T*dt;
exit_surv       = zeros(tot_years,1);
b_surv          = zeros(tot_years,1);
z_surv          = zeros(tot_years,1);
surv_t          = zeros(tot_years,1);

n_surv=n_sim;inde_surviv=ones(n_sim,1);
for tt=1:tot_years
    t0= (tt-1)/dt+1;
    t1= tt/dt;
    n_died        = sum(sum(i_exit(inde_surviv>0,t0:t1)));
    exit_surv(tt) = n_died/n_surv;
    b_surv(tt)    = sum(sum(b_hist(inde_surviv>0,t0:t1).*z_hist(inde_surviv>0,t0:t1)))/sum(sum(z_hist(inde_surviv>0,t0:t1)));
    z_surv(tt)    = mean(mean(z_hist(inde_surviv>0,t0:t1)));
    n_surv        = sum(1-i_exit(inde_surviv>0,t0));
    surv_t(tt)    = n_surv/n_sim;
    inde_surviv   = b_hist(:,t1);
end
z_all     = mean(z_hist(:,end));


exit_prob_yearly  = exit_surv+delta*(1-exit_surv);
